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ABSTRACT 

We use a set of large cosmological N-body simulations to study the internal structure 
of dark matter haloes which form in scale-free hierarchical clustering models (initial 
power spectra P(k) oc k n with n = 0,-1 and —2) in an £1 = 1 universe. We find that the 
radius r^g in a halo corresponding to a mean interior overdensity of 178 accurately 
delineates the quasi-static halo interior from the surrounding infalling material, in 
agreement with the simple spherical collapse model. The interior velocity dispersion 
correlates with mass, again in good agreement with the spherical collapse model. 
Interior to the virial radius rns, the spherically averaged density, circular velocity 
and velocity dispersion profiles arc well fit by a simple two-parameter analytical model 
proposed by Navarro et al. (1995a). This model has p oc r _1 at small radii, steepening 
to p oc r -3 at large radii, and fits our haloes to the resolution limit of the simulations. 
The two model parameters, scalelength and mass, are tightly correlated. Lower mass 
haloes are more centrally concentrated, and so have scalelengths which are a smaller 
fraction of their virial radius than those of their higher mass counterparts. This reflects 
the earlier formation times of low mass haloes. The haloes are moderately aspherical, 
with typical axial ratios 1 : 0.8 : 0.65 at their virial radii, becoming gradually more 
spherical towards their centres. The haloes are generically triaxial, but with a slight 
preference for prolate over oblate configurations, at least for n — — 1 and 0. These 
shapes are maintained by an anisotropic velocity dispersion tensor. The median value 
of the spin parameter is A « 0.04, with a weak trend for lower A at higher halo mass. We 
also investigate how the halo properties depend on the algorithm used to identify them 
in the simulations, using both friends-of-friends and spherical overdensity methods. We 
find that for groups selected at mean overdensities ~ 100 — 400 by either method, the 
properties are insensitive to how the haloes are selected, if the halo centre is taken as 
the position of the most bound particle. 
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1 INTRODUCTION 

The structure of dark matter haloes is of fundamental im- 
portance to understanding the formation and evolution of 
galaxies and galaxy clusters. In hierarchical clustering mod- 
els, galaxies form by gas cooling and condensing in dark 
matter haloes, and clusters form by the gravitational aggre- 
gation of individual galaxies and galaxy groups. This evolu- 
tion is driven by gravitational instability via accretion and 
merging. Many studies have been made of the properties of 
haloes formed by dissipationless gravitational collapse using 
N-body simulations, which are ideally suited to investigate 



this inherently non-linear problem. Continuing advances in 
computer technology and codes have enabled ever larger 
simulations to be performed, allowing halo structure to be 
examined in ever increasing detail. 

Pioneering N-body studies of the structure of haloes 
formed by hierarchical clustering in an expanding universe 
were carried out by Frenk et al. (1985) and Quinn, Salmon 
& Zurek (1986). These concentrated on halo density pro- 
files, and found approximately flat rotation curves for Q — 1 
cold dark matter (CDM) models, corresponding to density 
varying with radius roughly as p oc r~ 2 . Frenk et al. (1988) 
made a detailed study of haloes in CDM models in flat, 
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open and closed universes, using a P 3 M code with 32 3 par- 
ticles, calculating rotation curves, shapes and angular mo- 
menta of haloes. Efstathiou et al. (1988) made a similar 
study of haloes in self-similar clustering models with Q = 1 
and scale-free initial conditions, meaning initial power spec- 
tra of the form P(k) oc k n , with n — 1, 0, —1 and —2. These 
studies, and those of Quinn, Salmon & Zurek (1986) and 
Zurek, Quinn & Salmon (1988), which used a combination 
of PM and PP methods and comparable numbers of par- 
ticles, found that the slopes of halo density profiles varied 
both with spectral index n and with fi. 

More recently, there have been studies with higher reso- 
lution and/or more particles. Warren et al. (1992) simulated 
Q, — 1 scale-free models using a tree-code and ~ 10 6 par- 
ticles. They concentrated on the radial variation of shapes 
and the distribution of angular momentum in their better- 
resolved haloes. Crone, Evrard & Richstone (1994) used 64 3 
particle P 3 M simulations to study the density profiles of 
haloes in a variety of cosmological models for scale-free ini- 
tial conditions. They fit these density profiles to power-laws. 
However, they only examined the most massive haloes in 
each simulation. Navarro, Frenk & White (1995b) used a 
hierarchical tree-code to simulate the formation of a small 
number of haloes of very different masses in the CDM model 
at high resolution, again to study the density profiles. Du- 
binski & Carlberg (1991) also simulated individual haloes 
at high resolution, but with a much cruder treatment of the 
tidal effects of material outside the halo. 

In this paper, we analyse the structure of haloes formed 
in self-similar clustering models, with = 1 and spec- 
tral indices n = —2,-1,0. Scale free initial conditions are 
both conceptually simple, and, over limited mass ranges, 
provide useful approximations to physical models (such as 
CDM) whose power spectra have a slope which varies slowly 
with mass. For example, in the standard CDM model, the 
effective slope is n ~ — 2 on the scale of galaxy haloes 
(M ~ 1O 12 M0) and n ~ —1 on the scale of galaxy clus- 
ters (M ~ 1O 15 M0). Such models also have the great ad- 
vantage that their self-similar scaling properties can be used 
to distinguish physical effects from artifical features intro- 
duced by the limitations of the numerical simulations. Our 
simulations used the P 3 M code of Efstathiou et al. (1985), 
with 128 3 w 2 x 10 6 particles. Unlike Warren et al. (1992), 
we employ periodic boundary conditions. As in the previous 
studies, we pay particular attention to halo density profiles 
and rotation curves, but we also study the radial variation 
of other dynamical properties, the departures from spherical 
symmetry and the dependence of these properties on initial 
conditions. We also investigate how these properties depend 
on halo mass. 

We are interested in haloes as objects in approximate 
dynamical (or virial) equilibrium. An important issue not 
addressed in detail in previous studies is where one should 
draw the boundary of the virialized region of the halo, and 
what group-finding method works best for partitioning a 
simulation into virialized objects. Previous studies have each 
just used a single method of identifying the haloes (although 
several methods have been tried), and presented results for 
that one method, without investigating how the halo proper- 
ties might depend on the method used. In the present paper, 
we consider several possible criteria for deciding where the 
boundary of the virialized region is. We also compare results 



obtained using two different group-finding methods, friends- 
of-friends (Davis et al. 1985) and the spherical overdensity 
method of Lacey & Cole (1994), for a range of overdensities. 

The plan of the paper is as follows: In Section || we 
define virial mass and length scales for haloes, and compare 
three analytical halo models which we later use to fit our nu- 
merical results. Section ^ describes the N-body simulations, 
and how we identify the groups within them and construct 
radial profiles of halo properties. In Section |^, we examine 
the scalings of bulk properties of the haloes with mass, and 
make a critical comparison of the different group finding al- 
gorithms. In Section ^, we examine in detail the spherically 
averaged halo profiles, which we compare with the analyti- 
cal models of Section Section |^ discusses the angular mo- 
menta of haloes. Departures from spherical symmetry are 
investigated in Section |^ Our understanding of these re- 
sults and their implications are discussed in Section [], and 
we conclude in Section H. 



2 HALO MODELS 

The analytical model most often used to describe the for- 
mation of dark matter haloes is the idealised collapse of a 
uniform, spherically symmetric overdense region. For fi = 1, 
collapse to a singularity occurs when linear theory would 
predict an overdensity of <5 C = 3/20 (127r) 2//3 ~ 1.69. If one 
argues that realistic amounts of substructure present before 
the collapse will lead to violent relaxation, then the virial 
theorem and energy conservation imply that after virialisa- 
tion the object will have a density of 18-7T 2 ~ 178 times the 
background density p, if the final equilibrium state is also 
a uniform sphere. Real dark matter haloes are not uniform 
spheres, but we will refer to the radius rng around the halo 
centre within which the mean density is 178p as the "virial 
radius". This terminology will be justified by the results 
presented in Section ^, where we will see that the radius 
J"i78 approximately demarcates the inner regions of haloes 
at r < ri78 which are in approximate dynamical equlibrium 
from the outer regions at r > rns which are still infalling. 
The corresponding "virial mass" and circular velocity are 
given by 

M 178 = — 178prf 78 , (2.1) 
and 



= ( 



GM 17S 



V ri78 



1/2 



(2.2) 



These quantities prove to be useful for characterising the 
global properties of realistic dark matter haloes. 

Below we describe three equlibrium models for the in- 
ternal structure and kinematics of DM haloes, which we will 
later compare with the properties of the haloes formed in a 
set of large cosmological N-body simulations. Th ese m odels 
are the singular isothermal sp here, t he Hernquist (199C) and 
the Navarro, Frenk & White (Il995a|) (hereafter NFW) mod- 



els. The latter two models differ from the isothermal sphere 
in having shallower density gradients in the inner parts and 
steeper gradients in the outer parts. All three models are 
intended to represent only the regions interior to the virial 
radius ri78. The Hernquist model has finite total mass, while 
isothermal sphere and NFW models both have infinite total 
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mass. The Hernquist and NFW models can both be con- 
sidered, in some sense, as only minor modifications of the 
isothermal sphere, as interior to the virial radius, the den- 
sity p is approximately proportional to r~ 2 over the region 
containing most of the mass. The density, circular velocity 
and velocity dispersion profiles of these models are compared 
in Fig. 1 for particular values of the scale parameters. 



2.1 The Singular Isothermal Sphere 

A standard model for dark matter haloes, motivated by flat 
rotation curves, is the singular isothermal sphere, for which 
the density profile is 



p(r) oc 



(2.3) 



Defining s = r/r 17s as the radial distance in units of the 
virial radius, we can rewrite this as 



p 



178 
3 



(2.4) 



The corresponding expression for the mass within radius r, 
scaled to the value at the virial radius, is 



M(s) 
— — = s 

Ml78 

The circular velocity is constant with radius, 

Vc{s) _ 



(2.5) 



(2.6) 



and, assuming dynamical equilibrium and an isotropic ve- 
locity distribution, the 1-dimensional velocity dispersion is 
also constant, 



Vl78 



= V2, 



(2.7) 



In their relatively low resolution CDM simulations, Frenk et 
al. (1985,1988) and Quinn, Salmon & Zurek (1986) found 
halo circular velocities which were essentially constant with 
radius down to the gravitational softening radius, consistent 
with this simple model. This model, truncated at the virial 
radius, been widely used as a description of the dark matter 
haloes formed in hierarchical clustering models, e.g. in the 
modelling of gravitational lensing statistics by Narayan & 
White (1988), and in the galaxy formation models of Kauff- 
mann, White & Guiderdoni (1993) and Cole et al. (1994). 



2.2 The Hernquist Model 



Hernquist (1990) presented a completely analytical model 
which closely approximates the de Vaucouleurs R 1 ^ 4 law 
that has long been used to fit the surface brightness pro- 
files of elliptical galaxies. The density profile has the form 
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Figure 1. Comparison of the density, p, circular velocity, V c , 
and velocity dispersion, a\P , profiles for the singular isothermal 
sphere (solid) and the isotropic Hernquist (dashed) and NFW 
(dotted) models. In the NFW model the parameter ajvj = 0.2, as 
advocated to model the density profile of a CDM galaxy cluster by 
Navarro et al. (1995a), and ajj = 0.43 in the Hernquist model so 
that both have circular velocity profiles that peak at r ~ 0.43ri78. 
The "virial radius" rim is marked by the vertical dotted line. 
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Defining a dimensionless radius s = r/rn$ as before, and a 
dimensionless scale radius oh = bii/ r i78, the density, mass, 
potential t[> and circular velocity can be written in dimen- 
sionless form as 



p 



178 



2 o H (l 



an) 



3 s(s + a H ) 



(2-9) 



M(s) 

Ml78 

n 2 78 

Vc(s) 
Vl78 



(1 



\2 2 
■ OH] S 



(s + a H ) 2 
(l + a H ) 2 



(s + an) 



(1 



(s + a H 



«Hj 1/2 

s 



(2.10) 
(2.11) 
(2-12) 



© 0000 RAS, MNRAS 000, 000-000 



4 S. Cole and C. Lacey 



The parameter an sets the scale at which the slope of the 
density profiles changes from r _1 at small r, to r -4 at large 
r. The circular velocity profile peaks at r = aHri78 (see 
Fig. 1). The main attraction of this model is that it has 
an analytical distribution function which can be expressed 
in terms of elementary functions. It has also been used by 
Dubinski & Carlberg (1991) to fit the density profiles of 
galactic haloes formed in their N-body simulations. 

Assuming dynamical equlibrium, the radial velocity dis- 
persion, oy(r), of the model can be obtained by integrating 
the Jeans equation 



ld_ 

p dr 



(p^(r)) + 2(3- 



dr ' 



(2.13) 



where (3=1 — ag(r)/a^(r). For the case of isotropic or- 
bits, the radial and tangential velocity disper sions are equal, 
&e{r) = cr T (r). For this case Hernquist (1990) finds 
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2.3 The NFW Model 

Navarro et al. (1995a) proposed an alternative analytical 
density profile, which they found to be an excellent fit to 
the dark matter haloes of galaxy clusters formed in their 
CDM simulations. Navarro et al. (1995b) found this profile 
to apply over a large range of masses from those of dwarf 
galaxies to galaxy clusters. The density profile has the form 



p(r) oc 



r(r + 6n) 2 



(2.15) 



Defining on = ba/rns, the density, mass, potential and cir- 
cular velocity can be written as 
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where, for compactness, we have defined 

/(<IN) = (ln(l + l/a N )-l/(l + a N ))- 

Like the Hernquist profile, p oc r _1 as r — > 0, but the asymp- 
tote at large r is now r~ 3 rather than r~ 4 , which leads to 
a logarithmically divergent total mass. The parameter on, 
again, defines the scale of the transition between the in- 
ner and outer slopes, which occurs more gradually in this 
model than in the Hernquist model (see Fig. 1). The slope 
of the density profile is din p/ din r — —2 at r — a^ms- 
The circular velocity peaks at r — 2.16aNri78. Also shown in 
Fig. 1 is the 1-dimensional velocity dis persi on, a 1 . , obtained 
by numerically integrating equation ( (2.13[ ), again assuming 
isotropic orbits and dynamical equilibrium. 



3 SIMULATIONS AND GROUPS 



3.1 Simulations 



The simulations we have analysed are the three scale- 
free simulations which were used to study the statis tics o f 
mergers in an SI = 1 universe by Lacey & Cole (1994). 
The simulations were performed using the high resolution 
par ticle-p article-particle-mesh (P 3 M) code of Efstathiou et 
al. ( p.985p with 128 3 » 2x 10 6 particles. The long-range force 
was computed on a 256 3 mesh, while the softening parameter 
for the short-range force was chosen to be r\ = 0.2(L/256), 
where L is the size of the (periodic) computational box. The 
corresponding potential can be approximated by a Plum- 
mer law with softening e ss r//3 = L/3840. The initial power 
spectra of the simulations were power laws, P(k) oc k n , with 
n = —2,-1 and 0. More details are given in Lacey & Cole 



( p.994p . 

In this paper we analyse almost exclusively the last out- 
put from each of the three simulations. At this output the 
three simulations have expanded by a factor 27.8, 13.3 and 
6.3 and the characteristic mass is M» = 266, 447, and 46.8 
particles for n = 0, —1 and —2 respectively. Here we define 
the characteristic mass, M,, such that the r.m.s. linear den- 
sity fluctuation in spheres containing mass Af* is 1.69 - this 
being the linear theory density corresponding to the collapse 
to a singularity of a uniform spherically symmetric pertur- 
bation. To test for the effects of resolution we occasionally 
analyse a set of earlier outputs from these simulations at 
which the values of M» were a factor of 4 smaller. 



3.2 Group Finders 

We wish to study the intrinsic properties of non-linear self- 
bound structures that form via gravitational instability in 
the N-body simulations. For this reason we have identified 
these structures using different group finding algorithms to 
ensure that the properties we quantify are intrinsic to the 
structures and not artifacts of any particular group finding 
algorithm. Thus to identify the groups in the N-body simu- 
lation we employed two different group finding algorithms, 
the standard friends-of-friends (hereafter FOF) algorithm of 
Davis et al. (1985) and the spherical overd ensity (hereafter 
SO) algorithm described in Lacey & Cole (1994). 

FOF groups are constructed by linking together all pairs 
of particles whose separation is less than b times the mean 
inter-particle separation. This results in groups bounded by 
a surface of approximately constant density, p/p w 3/(27r6 3 ). 
This elegantly simple algorithm has been used extensively 
in previous analyses of N-body simulations, usually with b — 
0.2 (e.g. Frenk et al. 1988, Efstathiou et al. 1988). It succeeds 
in picking out most of the groups one can identify by eye, but 
occasionally can join together two or more distinct density 
centres that are linked by a tenuous bridge of particles. We 
will denote groups defined using this algorithm by FOF(6) . 

The SO algorithm first ranks all the particles in the 
simulation by their local density, computed from the dis- 
tance to their 10th nearest neighbour. Then, starting with 
the densest, a sphere is grown around this centre until the 
enclosed density drops below some threshold, up. The group 
centre is redefined to be the centre of mass of the selected 
region and the process of growing the sphere and selecting 



© 0000 RAS, MNRAS 000, 000-000 



The Structure of Haloes 5 



Table 1. The distribution of SO(178) groups by mass in the 
three simulations. The last line gives the mass M* in terms of the 
number of particles. 
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a new region repeated iteratively until the group centre and 
group membership converge. Particles assigned to a group in 
this way are removed from the density ranked list and not 
considered further. Then in the same manner each of the 
remaining ungrouped particles in the ranked list are used 
as group centres. Finally, small groups inside or overlapping 
with larger groups are merged into the larger object. This 
algorithm selects spherical regions whose average density is 
equal to up. A similar method, with a mean overdensity 
of 180, has been used by Warren et al. (1992). We denote 
groups found in this way by SO(k). 

3.3 Halo Profiles 

In constructing halo profiles we have chosen to include all 
the particles surrounding the group and not only those as- 
signed to the group by the group finding algorithm. In each 
of the groups defined above we locate the particle with the 
most negative gravitational potential energy. The potential 
is calculated using only the group particles, and assuming 
that the potential due to each particle is — l/(r 2 + e 2 ) 1//2 , 
with e = ri/3. We then use this as the centre of the group and 
inflate a sphere around this centre until the enclosed density 
drops below 178p. The radius of this sphere and the mass 
it encloses define ri78 and Mng for the group. We use ri78, 
M178 and the corresponding Vn$ to scale the halo profiles. 
These choices allow a fair comparison of groups of differ- 
ent masses and groups located using different algorithms. 
The question of which algorithm best divides the simula- 
tion into discrete virialized groups then becomes largely a 
matter of which returns a group mass close to the virialized 
mass, which we will see is close to our adopted definition of 
the virial mass Mng. The relation between Mng and the 
original group mass, A/ groU p is studied in Section ^| 

For each group centre we estimate a comprehensive se- 
lection of the properties of the surrounding material binned 
both in spherical shells and cumulatively within spheres. 
These properties include: the density, circular velocity, ra- 
dial velocity, radial and tangential velocity dispersions, ro- 
tation velocity, angular momentum, axial ratios of the mo- 
ment of inertia tensor and virial ratio. All internal veloci- 
ties are computed relative to the centre of mass velocity of 
the sphere of radius rns- The radial bins were chosen to 
be equally spaced for the innermost 10 bins and then with 



logarithmically increasing widths out to an outer radius of 
10ri78. The spacing of the inner bins was chosen such that 
the innermost bin, which in general is the least populated, 
typically contained 10 particles. When computing profiles 
averaged over all haloes within a mass range, the bin radii 
were scaled with ri78. We examine individual profiles con- 
structed in this way, and also mean profiles of the scaled 
quantities (p/p, V c /Vi78 etc.) averaged over groups within 
mass bins. The mass bins were chosen to cover a factor 2 in 
mass. With the exception of the density profile, the mean 
profiles evaluated in spherical shells are weighted means. We 
gave equal weight to all shells in the individual profiles that 
contained 10 or more particles and zero weight to those con- 
taining fewer than 10 particles. This procedure was adopted 
to avoid low particle numbers producing a large scatter in 
profiles such as the velocity dispersion as a function of ra- 
dius. For each of the mean profiles we also accumulated the 
group-to-group dispersion at each radius. These dispersions 
can be used to place errorbars on each of the profiles which 
reflect the group-to-group variation. 

Table m gives the distribution of masses in each of the 
three simulations for centres identified by SO(178). The 
large numbers of groups enable the mean profiles in each 
mass range to be determined very precisely. This then en- 
ables the dependence of halo profiles on mass and spectral 
index to be examined in some detail, but care must be taken 
to assess how resolution, which in terms of rm worsens with 
decreasing Mng/M*, affects these profiles. It is for this rea- 
son that we will also examine the profiles of groups extracted 
from an earlier output of the simulations when M* was a 
factor 4 smaller. 



4 BULK PROPERTIES 

We first look at the basic global properties of individual 
groups selected using the various group finders. As we will 
see in Section || when we study the dynamical structure of 
the haloes, the radius 7-178 corresponds well with the bound- 
ary at which there is a transition from general infall to ap- 
proximate dynamical equilibrium. Thus the virial mass of 
the halo is well approximated by M178. It is therefore in- 
teresting to see how the group mass defined by the various 
algorithms compares with Mi 78, as this acts as a measure 
of how well the particular group finders succeed in parti- 
tioning the simulation into discrete virialized systems. We 
see in Fig. ^ that for each group definition the group mass, 
Afg r0 up, correlates strongly with Mi 78, but the offset of the 
correlation and the scatter around it vary. The offset be- 
tween Mgroup and Afi78 is easily understood. Groups found 
using high k or small b are typically denser and less massive 
than those found using small k or large b, and so they are 
the central regions or sub-clumps of the less dense groups. 
Note that for FOF the offset is smallest for FOF(0.2). The 
scatter in each relation is of more interest. As is to be ex- 
pected from the similarity of their definitions, the scatter 
between Mg roup and M178 is smallest for SO(178) groups. In 
this case, the essential difference in the definitions of Mg roup 
and M178 is the choice of group centre. For SO(178) the 
group is centred by its centre of mass, while M178 is evalu- 
ated around the centre defined by the particle in the SO (178) 
group with the most negative gravitational potential energy, 
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Figure 2. Scatter plots of M gr0 up versus Mng for groups identified in the n = — 1 simulation using each of the group finders SO(400), 
SO(178), SO(IOO), FOF(0.15), FOF(0.2) and FOF(0.3). The mass M* equals 447 particles for this simulation. The straight indicates the 
locus Afg roU p = Mi7g. This line is omitted from the SO(178) plot where it would obscure most of the data points. The offset given in each 
panel is the median of log 10 (Afg roU p/Mi78). Also specified is the width of the distribution in \og 10 (M gToup / Mi7g) measured between the 
10th and 90th centiles. 



where the potential energy is evaluated using only the group 
particles. This difference is generally not significant as can 
be seen by the fact that many of the points fall exactly 
on the line M gTOUp = Mi 78. For SO groups identified with 
other values of the density contrast, ft, the scatter grows. 
For the dense SO (400) groups there is a tail of objects with 
M178 3> Mg roU p, for which the SO(400) group is simply the 
core of a larger virialized halo. For FOF(0.2) or (0.3) the 
scatter about the mean correlation is asymmetric with a 
tail of haloes with Mg roup > Mi 78. As the linking length b 
is increased this tail becomes more pronounced, extending 
to M g roup > 2Mi78 . These outliers arise when two or more 
distinct density centres become linked by tenuous bridges of 
particles. 

Fig. ^ shows the correlation of the 1-dimensional veloc- 
ity dispersion averaged within ri78 with the mass Mi 78. The 
correlation agrees well with the expectation of the isother- 



mal sphere model (Section 2.1), a 



Vi 78 /\/2 oc MjVg 3 , 
with an r.m.s. scatter about this mean relation of only 
10-15%. For the groups selected at high density contrast, 
particularly SO (400), a tail of groups appears with a™ > 
Vi.78/%/2. In these cases, it seems likely that the sphere of 
radius ri78 is actually only part of a significantly larger viri- 
alized structure. 

Groups selected at low density contrast, particularly 



those located by FOF, tend to be composities of distinct 
virialized structures, while, on the other hand, groups iden- 
tified at higher density contrast are often only small clumps 
within larger structures. Overall the FOF(0.2) and SO(178) 
group masses and M178 correlate well. Thus both group find- 
ing algorithms are good candidates for selecting virialized 
structures. The disadvantage of the SO(178) algorithm is 
that it forces spherical boundaries on the groups when we 
shall see that groups are typically triaxial in shape. The 
FOF(0.2) algorithm does not force any particular geometry 
on groups, but because it is only sensitive to the local den- 
sity it is prone to occasional merging of separate virialized 
systems that are linked by tenuous bridges. 



5 SPHERICALLY AVERAGED PROFILES 

We now examine the spherically averaged radial profiles of 
haloes and compare them to the analytical models described 
in Section ^| 



5.1 Massive Haloes 

Figs. ^| ^ and ^| show various mean profiles (averaged over 
mass bins) around the centres of the massive SO(178) and 
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Figure 3. Scatter plots of mean 1-dimensional velocity dispersion, cr^ D , within the virial radius, ri78, versus the mass Mi 78 for groups 
identified in the n = -1 simulation using each of the group finders SO(400), SO(178), SO(100), FOF(0.15), FOF(0.2) and FOF(0.3). The 
mass M, equals 447 particles for this simulation. V* is V178 for a halo with Mi 78 = M,. The straight indicates the locus cr^ D = Vng/2 1 / 2 , 
which is the prediction of the isothermal sphere model. 



FOF(0.2) groups from each of the three simulations. Here, 
and in all subsequent plots of radial profiles, two dotted ver- 
tical lines mark the force-softening scale r\ and virial radius 
ri78. Note that the effective softening for the P 3 M code is of- 
ten taken be e = rj/3. The upper two panels show the mean 
circular velocity, V c = (GM(< r)/i") 1 ^ 2 , and density pro- 
files. The density profiles flatten and the circular velocities 
turn up at large radii, r 3> rns, where material not physi- 
cally bound or associated with the group becomes included 
in the profile. Thus the density tends to the mean density 
at large r and the circular velocity becomes proportional to 
radius. Also shown on these figures are the corresponding 
Hernquist and NFW profiles which have been fitted to the 
SO(178) mean density profile over the range rj < r < ri78. 
The fits were performed by converting the group-to-group 
dispersion of the density profiles at each radius into an ef- 
fective error on the mean profile. The scale parameters cm 
and on of the Hernquist and NFW models were then varied 
to achieve minimum-^ 2 fits. Visually, both models are excel- 
lent fits to the halo density profiles between the virial radius, 
ri78, and the force-softening scale 77, for each value of n. The 
NFW model in particular continues to be a good match to 
the halo profiles at larger radii, where we will see that the 
halo is not yet in equilibrium, and at smaller radii where res- 
olution is likely to have affected the halo structure. A more 



telling comparison is that between the model and N-body 
circular velocity profiles. Note that by definition these are 
all constrained to pass through the same point at ri78. We 
see that the Hernquist model changes slope too abruptly to 
match that of the haloes. The NFW profile is a significantly 
better match. 

The lower panels of Figs. §, I and I show various mean 
kinematic properties of the haloes. These are the rotation 
speed Kot in spherical shells around the axis defined by the 
angular momentum vector of the shell, the mean radial ve- 
locity Vt, radial velocity dispersion <r r = {(V T — (14}) 2 ) 1//2 , 
and velocity dispersion anisotropy ag(r)/a T (r), where a 2 = 
{Vg} is the tangential velocity dispersion. In each case the 
velocities are defined with respect to the centre of mass mo- 
tion of all the material interior to radius ri78. To a first 
approximation the rotation velocity is small and the veloc- 
ity dispersion isotropic and almost independent of radius, 
all in agreement with the simple isotropic isothermal sphere 
model. A more detailed inspection shows that Vrot/ov ~ 0.3, 
which although appreciable, implies that the rotation is not 
a significant contributor to the dynamical support of the 
halo. Interior to ri78, cre/cr? is slightly less than unity, indi- 
cating a tendency for radial orbits, particularly at r ~ ri78. 
Outside the virial radius ag/a r increases and slightly over- 
shoots unity before settling to the isotropic value, ag /a T — 1, 
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Figure 4. The upper and middle panels show mean circu- 
lar velocity and mean density profiles of all groups with 2 < 
Afire /-W* < 4 (532 < N < 1064) in the n = simulation located 
by both SO(178) and FOF(0.2). In this case the lists of group cen- 
tres defined by SO(178) and FOF(0.2) were identical. The smooth 
curves in these two panels are the Hernquist (dashed) and NFW 
(dotted) models which have been fitted to the density profiles 
over the range between the vertical dotted lines which mark the 
force-softening scale rj and virial radius ri7g. The lower panel 
shows a variety of measures of the dynamical state of the haloes. 
The lowest solid curve is the mean radial velocity, V r , as a func- 
tion of radius. The next solid curve is the mean rotation speed, 
Vrot • The upper dot-dashed curve is erg /<r r , which measures the 
anisotropy of the velocity dispersion. The remaining solid curve 
traces the mean radial velocity dispersion as a function of ra- 
dius. The smooth dashed and dotted curves are respectively the 
Hernquist and NFW models for the velocity dispersion assuming 
isotropy, ag /<r r = 1 . 
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Figure 5. As Fig. ^, but for groups from the n = — 1 simulation 
with masses in the range 4 < Myjs/M* < 8 (1788 < N < 3576). 
Here the lists of group centres produced by the SO(178) and 
FOF(0.2) algorithms were slightly different and thus two slightly 
different curves are visible for each profile. 



at large radii. Overlaying the velocity dispersion profiles are 
the model predictions of the fitted Hernquist and NFW mod- 
els, computed assuming spherical sym metry and og /o> = 1 
(P = 0), from the Jeans equation (2.13). We see that, within 
the virial radius, both match well the slow curvature in the 
corresponding halo a T profile. 

The lowest curve in the lower panels of Figs. and ^ 
is the mean radial velocity, V t of material in shells. These 
curves are the main evidence we have that the radius rns 
is an accurate characterisation of the virial radius of the 
haloes. In each case Vt is approximately zero interior to 
ri78, consistent with a dynamically relaxed system. Between 
ri78 and 4rn$ V T < indicating material in the process of 
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Figure 7. The density and radial velocity profiles of a selection 
of 5 groups taken from the n = — 1 simulation with masses in the 
range 8 < M 178 /M* < 16 (3576 < N < 7152). 



falling onto the central halo. At larger radii, V T begins to 
increase passing through V r = at the turn-around radius 
and matching onto the general Hubble flow at large r. The 
location of the turn-around radius is in good agreement with 
the spherical collapse model which predicts a mean density 
interior to the turn-around radius of (37r/4) 2 w 5.55. Simi- 
lar results for the radial velocity were found by Crone et al. 
(1994). 

The typical group-to-group variation of the radial ve- 
locity pattern and the corresponding density profiles is il- 
lustrated in Fig. which shows the profiles for 5 groups 
taken from the n — — 1 simulation with masses in the range 
8 < M/M* < 16. Although the profiles can be significantly 
perturbed by substructure in the individual haloes, the tran- 
sition from infall to quasi-static equilibrium at r w ri78 is 



quite clear in the majority of cases. But the radial velocity 
profiles are too noisy to make them useful for determining 
the boundary of the virialized region for individual haloes. 

Another measure of dynamical equilibrium within 
haloes is provided by the virial ratio 2T/\W\. The expec- 
tation for an isolated system in dynamical equilibrium is 
for the total binding energy to equal twice the total kinetic 
energy, 2T = \W\. As an example, Fig. ^ plots the mean 
value of the virial ratio 2T /\W\ evaluated within spheres for 
groups from the n = — 1 simulation with masses in the range 
8 < Mns/M* < 16. Here T is the kinetic energy of the mate- 
rial within a sphere of radius r and W the self-gravitational 
binding energy of the same material, i.e. neglecting the the 
effect on the gravitational potential of all material at larger 
radii. The ratio 2T/\W\ should approach unity at the bound- 
ary of the virialized region if surface terms in the virial theo- 



© 0000 RAS, MNRAS 000, 000-000 



10 S. Cole and C. Lacey 



V V 


i 
1 


1 1 1 1 1 1 


1 1 1 


1 1 1 1 


" \ \ 






n=-l 




; \ 














8<M 17fl/ 


7 M t <i6: 




\ 

. \\ 








- 


'• V 
'•-\\ 






- 






















i 


, 1 , 


1 1 1 1 1 1 




II 



2.5 



-1.5 



-0.5 
lo gi( 







0.5 



1 



r/r, 



10\ L I - 1 178/ 

Figure 8. The solid curves shows the mean virial ratio, 2T/\W\, 
evaluated within spheres, as a function of radius, for haloes with 
masses in the range 8 < M 17s /M* < 16 (3576 < N < 7152) from 
the n = — 1 simulation. The dashed and dotted curves show re- 
spectively the corresponding profiles for the Hernquist and NFW 
models that were fitted to the mean density profiles of these 
groups. 



rem vanish. We find 2T/\W\ approaches unity only slowly at 
large radii and is significantly greater than unity at r = rng. 
However this same profile is tracked quite well by the same 
quantity evaluated from the fitted analytical Hernquist and 
NFW models, which when integrated over all the mass ac- 
curately obey the virial theorem. Thus, although the haloes 
are believed to be in reasonably good dynamical equilibrium 
within r < rn$, the ratio 2T/|W| is not useful for defining 
the boundary of the virialized region. 



5.2 Mass dependent Trends 

In the previous section we saw that the NFW model pro- 
vides an excellent description of the spherically averaged 
structure of the more massive haloes formed in each of the 
three simulations. We now turn to how the halo structure 
depends on mass. The density and circular velocity profiles 
for a range of masses in each of the three simulations are 
shown in Figs. [] and [lo| Also shown are Hernquist and NFW 
models which have been fitted to the density profiles over 
the range rj < r < ms- Over this range, the NFW model 
is an excellent description of the halo density profiles for 
almost the full range of masses investigated in each of the 
simulations. Very close to the softening radius rj the den- 
sity profiles of the most massive halos are slightly steeper 
than those of the NFW model. For lower masses the fit- 
ted NFW models also match the inner portion of the halo 
density profiles at r < 77, in all cases except the low mass 
groups of the n = —2 simulation. Although not apparent 
to the eye, the Hernquist models are formally worse fits to 
the halo density profiles than the NFW models. They tend 
to be marginally too steep in the outer parts of the haloes 
and too shallow in the inner regions. Nevertheless they are a 
considerable improvement over a simple isothermal model, 
which has p oc r~ 2 at all radii. 




1 

log(r/r 178 ) 



Figure 12. A comparison of the mean density profiles of SO(178) 
groups of a variety of masses in the n = — 1 simulation at two dif- 
ferent output times. The solid curves are the mean density profiles 
from the final output, and are the same as those plotted in the 
central column of Fig except here we plot log((r/ri7g) 2 p/(p)) 
so as to expand the vertical scale. The dashed curves are the cor- 
responding profiles from an earlier output in which the groups 
are less well resolved. At this earlier time M* = 112 particles, 
a factor four less than at the final time. The inner two vertical 
dotted lines mark the force softening scale r\ in the two cases. 



The differences between the NFW and Hernquist mod- 
els are more pronounced in plots of the circular velocity pro- 
files, Fig. Here we again see that the NFW models match 
well the profiles of the haloes over the entire region of the 
fit. In contrast the Hernquist models peak too sharply and 
underestimate the circular velocity in the central regions of 
the haloes. 

We now turn to the question of systematic variations 
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Figure 9. The mean density profiles of SO(178) groups of a variety of masses from each of the three self-similar simulations. The 
two dotted vertical lines mark the force-softening scale, r), and virial radius, ri78- The smooth curves are NFW (dotted) and Hernquist 
(dashed) model fits to the profiles in the range r\ < r < rijs- Note that M» equals 266, 447 and 46.8 in the n = 0, — 1 and —2 simulations 
respectively. 
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Figure 10. This figure shows the circular velocity profiles and model fits corresponding to Figure 



of the halo density profiles with mass, Mng/M,. Here we 
have to be very careful to distinguish between real physi- 
cal effects and the consequences of the limited resolution of 
the N-body simulations. If the finite number of particles or 
the force resolution were artificially supressing the density 
in the centre of the simulated haloes then we would expect 



the radius of the affected region to be directly proportional 
to r\. In contrast we see in Fig. |l^ that the peak in the circu- 
lar velocity profiles moves to smaller r/rirs for lower mass 
haloes while at the same time rj/rns, marked by the inner 
vertical line, increases. This is the first piece of evidence 
which gives us confidence that the form of the density and 
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Figure 11. This figure shows the mean profiles of kinematical properties corresponding to Figure |9j. The lowest solid curve in each panel 
is the measured radial velocity Vr, followed by rotation speed Vrot and radial velocity dispersion oy, and the upper dot-dash curve is the 
velocity dispersion anisotropy. The smooth dashed and dotted curves are the predictions of the isotropic Hernquist and NFW models 
respectively, based on the fits to the density profiles shown in Figure Bl 
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Figure 13. The variation of the scale parameters ajq and ajj of 
the NFW and Hernquist model fits to the halo density profiles. 
The solid symbols are the parameter values from the fits to the 
mean density profiles of the groups from the final output of each 
simulation. The circles, squares and triangles are for n = —2, —1 
and —2 respectively. M* for these outputs corresponds to 266, 477 
and 46.8 particles. The open symbols, plotted displaced slightly to 
the right, are the corresponding parameter values obtained from 
earlier simulation outputs when the characteristics masses, M,, 
were a factor 4 smaller. The error bars on each point indicate the 
±5cr range, where the formal la errors were obtained from the 
Ax 2 = 1 points of the x 2 -fits to the mean density profiles. 



circular velocity profiles is not predominately due to limited 
resolution. A more stringent test comes from analysing the 
haloes with the same values of Mirs/M* , but extracted from 
an earlier output of each simulation when M* was smaller 
by a factor 4, and therefore rj/rns larger by a factor 1.59. 
An example of this comparison is shown in Fig ^| which 
compares the mean density profiles of groups in the n = — 1 
simulation at the two output times. We note that the two 
density profiles agree quite accurately for radii greater than 
r\, the force softening scale of the less well resolved output. 
Fig jl3| shows the dependence of the scale parameters hn and 
bh on Mi 78 /M* for each of the three simulations and for the 
two output times. The fits were made to the binned mean 
density profiles over the range r\ < r < rirg using estimates 
of the error on the mean profile obtained from the group- 
to-group variance about the mean divided by the number 
of groups. Standard x 2 -fits were made, treating each data 
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Figure 14. The scatter in the virial ratio evaluated interior to 
ri7s versus mass for the n = — 1 simulation. The curves mark the 
locus of the 20th, 50th and 80th centiles of the distribution. The 
mass M* equals 447 particles. 



point as independent and formal la errors on the fit parame- 
ters estimated from the points adjacent to the best fit where 
Ax 2 = 1. The error bars on Fig [L3 show the ±5a bounds. If 
numerical resolution were the sole reason for the existence 
of the detectable break in the halo density profiles, then we 
would expect the scale lengths an or an to scale as Af _1//3 
at a given output time, and we would expect the values from 
the earlier output to be 60% larger than those from the final 
output. In fact they are typically less than 20% larger than 
their better resolved counterparts from the final output of 
each simulation. Thus, although the measured scale lengths 
must have some residual dependence on the numerical reso- 
lution, the detected trends with mass and spectral index are 
real. The scale lengths, on or an, increase with increasing 
mass M178/M*, and the steepness of this relation depends 
on the spectral index n. 

We also examined the mean profiles of velocity disper- 
sion, anisotropy, rotation velocity and radial velocity for 
haloes in all mass ranges. The results, shown in Figure [III 
are very similar for the different masses. In particular, the 
radial velocity always starts to go significantly negative (in- 
dicating the edge of the virialized region) at r ~ ri7g. 

The distribution with mass of the ratio (2T/\W\)ns, 
evaluated in spheres of radius ri78, is shown in Fig. [L4j, for 
the n = — 1 simulation. There is no tendency for the lower 
mass, less well resolved haloes, to be further from virial equi- 
librium than the larger mass haloes. In fact, the median 
value of (2T/| W|)i78 is 1.1 at low masses and increases grad- 
ually to approximately 1.2 at Mwg/M, = 10. This trend is 
completely consistent with the variation of the halo density 
profiles reported above. The ratio 2T/\W\ at the virial ra- 
dius can be computed for the NFW model. For cin = 0.1, 
which fits the density profiles of haloes of mass Mi 78 ~ M* , 
(2T/|W[)i78 = 1.13, while for hn = 0.2 which corresponds 
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Figure 15. The distribution of the spin parameter A178, eval- 
uated for material within a sphere of radius rng, against mass. 
The lines mark the locus of the 20th, 50th and 80th centiles of 
the distribution. Note that Af* equals 266, 447 and 46.8 in the 
n = 0, —1 and —2 simulations respectively. 



to MiTg ~ 10M*, (2T/|W|)i78 = 1-23. At all masses there 
is a tail of objects with larger values of (2r/|W|)i7g. These 
are probably examples of ongoing mergers such as depicted 
in Fig. 
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Figure 16. The alignment of the angular momentum within 
haloes. The angle 9j is the angle between the total halo angular 
momentum, computed from all the material interior to the virial 
radius ri78, and the angular momentum of material in shells or 
spheres. The curves give the mean value of cos(0j) for shells (solid 
curves) and spheres (dashed curves) where the averaging is done 
over all the haloes in the specified mass ranges. 



6 ANGULAR MOMENTUM 

The dimensionless spin parameter 

J\E\ 1/2 

where J, E and M are the total angular momentum, energy 
and mass is an important dynamical parameter that has 
been much studied in numerical simulations ( e.g. Barnes & 
Efstathiou 1987). Our results are in broad agreement with 



earlier work and have a median A ~ 0.04 at M = M*. Fig. |l5| 
shows, as a function of mass, the distributions of Ai78 eval- 
uated for spheres of radius ri78. The median, 20th and 80th 
centiles are indicated by the solid lines. They reveal no ob- 
vious variation in the distributions with spectral index, but 
a weak trend towards decreasing values of A178 with increas- 
ing mass. The earlier outputs from the simulations confirm 
that this trend is real, and not simply a result of the ef- 
fective numerical resolution varying with mass. The overall 
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distributions and trends with mass are in good agreement 
with those of Barnes & Efstathiou (1987), who studied the 
distribution of A for CDM and n = models, and those re- 
ported by Efstathiou et al. (1988), who investigated other 
scale free models. Efstathiou et al. also found a weak de- 
pendence of the median value of A on spectral index, but 
for a fixed Mi 78 /M* we find any such dependence to be ex- 
tremely weak. For example, for 1/2 < Mvrs/M* < 2, we find 
median values for A178 of 0.045, 0.043 and 0.044 for n = 0, 
— 1 and —2 respectively, with uncertainties of about 0.001. 
For 2 < M 178 /M* < 4, we find medians 0.034, 0.037 and 
0.038 with uncertainties of around 0.002. Efstathiou et al. 
may have found a stronger trend because by averaging over 
all the haloes in each simulation, they gave different weights 
to different ranges of M/M* for different n. 

The degree of alignment of the angular momentum 
throughout each halo is studied in Fig. hfl. For each halo 
we measured the angle 9j between the angular momentum 
vector of material in shells or spheres of radius r and the 
total angular momentum vector of the material within the 
virial radius rns ■ As a function of radius, Fig. |l^ shows the 
mean value of cos(8j) for shells (solid curves) and spheres 
(dashed curves) for haloes selected in various mass ranges 
from each of the three simulations. If the angular momen- 
tum were randomly orientated then (cos(0j)) = 0. We find 
some degree of alignment through the virialized region of 
the haloes, with a sharp drop in the alignment beyond rng. 
The estimate of the degree of alignment is sensitive to the 
noise in estimates of J. Thus, the true alignment may be 
somewhat better than indicated by this analysis. A more 
detailed study of the angular momentum distribution can 
be found in Warren et al. (1992). 



7 ASYMMETRY AND SUBSTRUCTURE 

In the previous analysis, by spherically averaging we have 
treated the haloes as if they were both smooth and intrin- 
sically spherical. In reality, the haloes are non-spherical and 
contain significant sub-structure. Figs. |l7ja-d show exam- 
ples of four groups selected from the n = — 1 simulation to 
illustrate the variety of morphologies. Each group is shown 
projected along each of the three principal axes of the in- 
ertia tensor of the material within the radius rrra, whose 
orientation is indicated by the dotted ellipse in each panel. 
The particles plotted are those identified as group members 
by the SO(178) algorithm. The contours show the projected 
surface density of all particles within 2ri78 of the group cen- 
tre. For these groups the outermost contours agree extremely 
well with the projected boundaries of the same groups when 
identified using the FOF(0.2) algorithm. 

The first two groups, Figs. |l7|a-b, are typical of the ma- 
jority of groups. There is a good correspondence between 
the FOF(0.2) and SO(178) group definitions and both have 
centres of mass that coincide quite well with the potential 
centre. The offset between these two centres is indicated in 
the upper panel of each figure. The remaining two groups 
are examples where the FOF(0.2) algorithm links together 
two distinct haloes which are in the process of merging. 
In such cases the offset between the potential centre and 
group centre of mass can be large for the FOF(0.2) group, 
but remains small for the SO(178) definition. Thus for the 
FOF(0.2) groups, there is a correlation between the mor- 
phology and the offset between the centre of mass and the 
potential centre. Approximately 20% of FOF(0.2) groups 
have doff > 0.3ri78. 

To quantify the distribution of halo shapes and their 
radial dependence within individual groups, we have studied 
the first and second moments of the mass distribution within 
spheres. 

The first moments define a dipole vector, whose first 
component is 

L x {r) = ^ z/iVsphere, (7.1) 
sphere 

where the summation is over all Adhere particles interior to 
radius r. The statistic that we study is a normalised mag- 
nitude of this vector, D(r) = |L(r)|/(r) r , where (r) r is the 
mass weighted mean value of r within the sphere. This def- 
inition implies that an isotropic distribution has D(r) — 0, 
while the maximum value of D(r) is unity corresponding to 
all the mass concentrated in one direction. 

The second moments of the mass distribution define the 
components of the inertia tensor. For example 

Ixy(r) = xy/N sphcrc , (7.2) 

sphere 

where the summation is again over all particles interior to 
radius r. If we denote the ordered eigenvalues of this tensor 
as a 2 , b 2 and c 2 , then the mass distribution can be charac- 
terised by an ellipsoid with axial ratios a : b : c. We study 
these axial ratios as a function of radius r, and also a nor- 
malised quadrupole statistic which may be written in terms 
of these eigenvalues as 
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a) 



b) 



Figure 17. A selection of 4 groups from the n = —1 simulation with masses in the range 8 < Mijg/M* < 16 (3576 < N < 7152), 
showing the variety of group morphologies and comparing the FOF(0.2) group definition with that of SO(178). The particles plotted are 
those defined as a group by the SO(178) algorithm. Each group is shown projected along the three principal axes defined by the inertia 
tensor of the material within the sphere of radius ms- In each panel the dashed circle marks the projected boundary of this sphere. The 
flattening and orientation of the inertia tensor is indicated by the dotted ellipses. The contours show the projected density of all particles 
within 2ri78. The contour levels are 0.5, 2 and 8 times the mean projected surface density within the rng circle. The outer connected 
contour corresponds quite accurately to the boundary of the group defined by the FOF(0.2) algorithm. The values of d ff given in the 
upper panel of each figure are the offsets of the potential centre from the group centre of mass for the SO(178) and FOF(0.2) group 
definitions. 
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Figure |l^. 



c) 

continued. 



.1/2 



(2(a 2 + b 2 + c 2 ) 2 - 6(a 2 fe 2 + b 2 c 2 + cV)V 
Q(r) = ^ '—^ (7 .3) 

where (r 2 ) r is the mass weighted mean value of r 2 within 
the sphere. An isotropic distribution has Q(r) — while the 
maximum value is Q(r) — v2- 

Figure |l8| shows the mean axial ratios and quadrupole 
and dipole statistics for a selection of well resolved haloes 
from each simulation. For fixed M\-re,/M* poorly resolved 



d) 



haloes from the earlier output time of our simulations have 
noisier and larger quadrupoles. For all mass ranges the halo 
shapes become increasingly aspherical - small values of b/a 
and c/a and large values of Q(r) and D(r) - at radii smaller 
than the force softening scale r\. At larger radii, for haloes 
containing more than approximately 100 particles within 
ri78, we find that these profiles have very little systematic 
dependence on mass. One can also see from Fig. [L8| that 
the profiles also have negligible dependence on the spectral 
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Figure 18. For each of the three simulations the curves in the upper panels show the axial ratios b/a (solid) and c/a (dashed) of the 
inertia tensor of material in spheres. The curves shown for n = are for masses 1 < Mvrs/M* < 2 and 4 < Mns/M* < 8. For n = —1 
curves for the two mass ranges 2 < M\ig/M t < 4 and 8 < M\i&/M* < 16 are shown. For n = —2 the mass ranges are 8 < Mus/M* < 16 
and 32 < Mm/M r < 64. The curves in the lower panels show the corresponding dipole, D(r) (dashed), and quadrupole Q(r) (solid) 
values. As before, the vertical dotted lines mark the force-softening scale r\ for the two mass ranges and the virial radius r\7g. Note that 
M* equals 266, 447 and 46.8 in the n = 0, —1 and —2 simulations respectively. 



slope of the initial conditions. The typical mean axial ratios 
a : b : c are 1 : 0.8 : 0.65 at the virial radius fire- The ax- 
ial ratio profiles indicate that the haloes gradually become 
more spherical with decreasing radius until one approaches 
the force softening scale. The dipole and quadrupole statis- 
tics are slowly increasing throughout the range r\ < r < rm, 
with typical values of 0.15 and 0.25 at the virial radius. Be- 
yond the virial radius, the random distribution of the mate- 
rial falling towards the halo results in a sharp decline in the 
axial ratios b/a and c/a and a corresponding sharp increase 
in both D(r) and Q(r). 

The distribution of axial ratios a : b : c at the virial 
radius is shown in Fig.^. This figure combines all haloes 
with masses Mn$ > M*, but splitting the samples reveals 
no significant mass dependence. For each spectral index the 
distribution of axial ratios is broad, with all haloes being 
to some degree triaxial. For spectral indices n = — 1 and 
prolate configurations (c/b < b/a) axe mildly favoured over 



oblate (c/b > b/a), but no such trend is evident in the case 
of n = -2. 

Our conclusions about axial ratios are in general agree- 
ment with those of Efstathiou et al. (1988). Our values for 
the typical axial ratios also agree with Warren et al. (1992), 
when for the latter we compare to the values they calculate 
from the inertia tensor within spheres. Warren et al. also 
use a method of fitting ellipsoids to the density distribution, 
from which they find a trend for haloes to become more 
spherical with increasing radius, opposite to the trend we 
find. 
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Figure 19. The distribution of the axial ratios b/a and c/b for the 
interia tensor of the mass interior to rirs for SO(f 78) groups with 
masses Mng > M* from each of the three simulations. The locii 
of prolate (c/b = 1) and oblate (b/a = 1) shapes are indicated. 



8 DISCUSSION 

To the limit of the resolution of our simulations, the densities 
of dark matter haloes diverge as one approaches their cen- 
tres. As a first approximation, the density profiles between 
the gravitational softening scale and the virial radius are fit 
by power laws, in agreement with Crone et al. (1994). For 
Q — 1 and n = —2, the power law is close to p oc r -2 , con- 
sistent with flat rotation curves, but the effective power-law 
slope steepens significantly with increasing n and with de- 



creasing mass. However, the density profiles have significant 
curvature in the logp - logr plane, steepening with increas- 
ing radius. Better 2-parameter fits than a power law p oc r~ a 
to the spherically averaged density a nd cir cular velocity pro- 
files are provided by the Hernquist (199C) or Navarro et al. 
(1995a,b) analytical models, which each have a scale radius 
b , and thus a logarithmic slope which changes with radius. 
Both of these models have density profiles which are shal- 
lower than r~ 2 at r <C b and steeper at r 3> b. The Hern- 
quist model changes slope too abruptly to match the sim- 
ulated haloes over the full range of resolved scales, but the 
NFW model is an excellent fit to all our mean halo profiles. 
The NFW model also reproduces other dynamical proper- 
ties of the simulated haloes interior to their virial radii rng. 
Ignoring the small radial anisotropy of the halo velocity dis- 
persion tensor, the slow variation of velocity dispersion with 
radius is well matched by that in the isotropic NFW model. 
Once the radius ri78 has been measured for a halo, the NFW 
model has only one free parameter, a dimensionless scale ra- 
dius on = b^/rirs- It is remarkable that such a simple one- 
parameter function is able to fit the density profiles of haloes 
over a wide range of masses and initial conditions. Formally, 
the Hernquist and NFW models imply density varying as 
p oc r" 1 at very small radii (r <C b), but even in our best re- 
solved haloes, our profiles only extend down to t/6n ~ 0.3, 
where the slope for the NFW model is approximately r -1 ' 5 , 
so we have not proved that the asymptotic r _1 dependence 
is correct. 

The dependence of the dimensionless scale radius, a, on 
halo mass and on the slope of the initial power spectrum 
is of particular interest. Lacey & Cole (1993,1994) demon- 
strated that higher mass haloes have higher accretion rates 
and more recent formation epochs than their lower mass 
counterparts, where the formation time is defined as the 
time at which half the present mass of a given halo had al- 
ready been assembled. Also, for a given Mm/M*. the halo 
formation time is earlier for increasing spectral index n. If 
the inner regions or "core" of the final halo formed mostly 
from the material which assembled earlier, then they should 
reflect the background density of the universe at the time 
they were assembled, and be denser for haloes which formed 
earlier. Therefore one predicts denser, more concentrated, 
cores for haloes with earlier formation times. This is pre- 
cisely the trend we find in the simulated haloes. For fixed n, 
the lower mass haloes are more concentrated, with smaller 
values of a, while at fixed mass, a decreases with increasing 
n. Navarro et al. (1995b) investigate the first of these trends 
for the CDM model, and find a quite good proportionality 
between the core density and the density of the universe at 
the formation time. 

The dependence of the density profiles and rotation 
curves on the spectral index is also in qualitative agreement 
with the results of Quinn et al. (1986), Efstathiou et al. 
(1988) and Crone et al. (1994), who found steeper profiles 
for power spectra with more small scale power, i.e. increas- 
ing n. Crone et al. (1994) found no significant dependence 
of the density profiles on halo mass. However, they exam- 
ined a smaller range of halo masses than in the present pa- 
per. They also applied a correction procedure intended to 
compensate for the effect of the finite force resolution in 
their simulations. This correction made all their haloes have 
steeper core density profiles. We tested this procedure on 
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our haloes with similar results. However, we did not find 
that it improved the agreement between the halo profiles 
with the same Mirs/M* from two different output times. 
Thus, although the procedure alters the density profiles, it 
is not clear that it is usefully correcting for all the effects of 
limited resolution, which include not only force resolution 
but also finite particle number and limited range of scales 
in the initial conditions. Thus we preferred to leave the den- 
sity profiles uncorrected and instead fit the analytical mod- 
els only over the range r) < r < rn$, where comparison of 
the different output times showed the density profiles to be 
robust. 

We can also compare our results with the analytical 
model of secondary infall around density peaks developed 
from the work of Gunn & Gott (1972) by Hoffman & Shaham 
(1985). Hoffman & Shaham predicted power- law halo den- 
sity profiles, with logarithmic slope equal to — 2 for n < — 1 
and — (9 + 3n)/(4 + n) for n > — 1. We find that the profiles 
depart from power laws, though the average profile slopes 
are similar to the Hoffman-Shaham values over the regions 
where we can measure them, and Hoffman-Shaham predicts 
the correct trend of steeper profiles with increasing n. How- 
ever, the Hoffman-Shaham model also predicts that the pro- 
file shapes should be independent of halo mass, which dis- 
agrees with our N-body results. 

We have investigated how best to draw the bound- 
ary between the region of a halo which is virialized and 
in approximate dynamical equlibrium, and the outer parts 
which are still infalling. This is particularly important if one 
wants to define a total virialized mass for the halo, to use 
in comparing with theoretical models which express results 
in terms of a total halo mass. An example of the latter is 
the Press-Schechter model for the mass function of haloes 
(Press & Schechter 1974), and its extension to statistics of 
halo mergers (Lacey & Cole 1993,1994). Our starting point 
was the spherical collapse model, which for the collapse of 
a uniform sphere in an Q = 1 universe, predicts that the 
virialized region lies interior to a radius ri7g for which the 
mean density is 178 times the background value. We then 
find that the mean radial velocity profiles plotted in the di- 
mensionless form V r /Vns vs. r/ri7s all look very similar for 
different n and different Mns/M*, and all show a break at 
r « rn S from a quasi-static interior with V r ~ to an in- 
falling exterior (for 1 < r/rvrs < 4) with V T < 0. Thus, 
the radius rns , and the corresponding enclosed mass Mi 78 , 
do indeed provide good characterizations of the radius and 
mass of the virialized region. The shape of the mass distri- 
bution also becomes significantly more non-spherical beyond 
ri78, consistent with a break between dynamically relaxed 
inner regions and unrelaxed outer regions, and supporting 
the choice of ri78 as the virial radius. On the other hand, the 
virial ratio 2T/\W\ in spheres varies smoothly through ri78, 
and at ri78 exceeds the value 2T/\W\ = 1 for an isolated 
object in dynamical equilibrium by 10-20%, principally be- 
cause of the confining effect of the external infalling material. 
Thus, the value of 2T/\W\ does not provide a good criterion 
for deciding where the virialized region ends. 

Crone et al. (1994) followed a similar procedure to us, 
but instead chose an overdensity of 300 to define the edge of 
the virialized region. In practical terms, the difference from 
our choice is not that large. For the haloes in our simulations, 
the dimensionless scale radius spans the range 0.05 <. on 



0.3. For this range, the NFW model gives r3oo/ri78 ~ 0.8 
and M300/M178 « 0.8 — 0.9, where r3oo and M300 are defined 
analogously to ri7 8 and Mns- 

The halo properties we analysed are largely independent 
of the algorithm used to identify the groups, for the SO 
method with mean overdensity in the range k = 100 — 400, 
and for the FOF method with linking length in the range 
b — 0.15 — 0.3, corresponding to a local overdensity in the 
range 20 — 140. The halo profiles are insensitive to exactly 
which particles are identified as group members because we 
use the minimum in the gravitational potential to define the 
halo centre, and include all the surrounding material when 
constructing the halo profiles. 

What group-finder partitions the particle distribution 
into objects that best correspond to virialized haloes accord- 
ing to the results above? Not surprisingly, the raw groups 
found by the SO (178) algorithm mostly match very well with 
the material interior to ri78 around the halo centre, defined 
as the minimum in the gravitational potential. The differ- 
ence is that the raw SO(178) groups are defined so that 
the centre of mass of the group coincides with the centre of 
the sphere. The FOF(0.2) groups mostly also match fairly 
well to the material within ri78 of the halo centre, as de- 
fined above, but occasionally the FOF(0.2) algorithm links 
two distinct haloes which are soon to merge, and in such 
cases the offset between the FOF(0.2) group centre of mass 
and the potential centre can be large. The masses of the 
FOF(0.2) groups mostly agree well with the virial masses 
Mi 78, with a spread around this of only about 20%, but 
with larger deviations for haloes in the process of merging. 
This provides some post hoc justification for the extensive 
use of the FOF(0.2) algorithm in computing halo properties 
and mass functions. 



9 CONCLUSIONS 

We have studied the structure of haloes that are formed 
in simulations of self-similar clustering models with n = 
0, — 1, —2 and Q = 1, over a wide range in halo mass. 

(1) We find that the radius ri78 about the halo cen- 
tre within which the mean overdensity is 178 accurately de- 
marcates the virialized interior of the halo, which is in ap- 
proximate dynamical equlibrium, from the exterior, where 
material is still falling in. The characterization of ri7s as 
the "virial radius" agrees with the prediction of the simple 
spherical collapse model. 

(2) For r < rns, the spherically averaged density, cir- 
cular velocity and velocity dispersion profiles are very well 
fit by the analytical model of Navarro et al. (1995a,b) with 
an isotropic velocity dispersion. This model has p oc r^ 1 
at r/rns <SC on, rolling over to p oc r~ a at r/rns ^ on, 
where on is a dimensionless scale radius. It should noted 
that even for out best resolved haloes the NFW model has 
p oc r~ at the softening radius, and thus we do not probe 
the asymptotic prediction of p oc r _1 . 

(3) The value of the dimensionless scale radius on cor- 
relates strongly with halo mass and with spectral index 
n. Haloes with mass M much less than the characteris- 
tic mass M* typically form much earlier than haloes with 
M 3> M*. This results in them being more centrally con- 
centrated (smaller on), with cores whose densities reflect the 
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higher mean density of the universe at their formation time. 
For a given M/M», haloes form earlier in models with more 
small-scale power (larger n), and also have denser cores. 

Despite the successes of the spherically symmetric 
model, haloes are not spherical. In these Q — 1 simulations 
most haloes show evidence of substructure and many arc 
clearly about to undergo a merger. 

(4) Haloes are generically triaxial, but with a slight pref- 
erence for prolate configurations, at least for n — — 1 and 0. 
Their inertia tensors (measured in spheres) have typical ax- 
ial ratios 1 : 0.8 : 0.65 at the virial radius ri78, becoming 
gradually more spherical towards their centres. There is no 
significant trend with mass or spectral index. Measurements 
of the dipole moments of the mass distribution around halo 
centres show them to have mean dipoles of order 10%, and 
similarly large quadrupole distortions. 

(5) Haloes are slowly rotating, with V ro t/cr w 0.3. The 
rotation velocity V ro t measured in spherical shells is approx- 
imately constant with radius, but the angular momentum is 
not well aligned between the central and outer regions. The 
median value of the spin parameter (measured interior to 
ri7 8 ) is A w 0.04 with a weak trend for lower A at higher 
halo mass. There is no significant trend of A with spectral 
index. 

(6) The properties we find are insensitive to the details 
of the group-finding algorithm used to select the haloes in 
the simulations, since once we have identified a group, we de- 
fine the halo centre from the minimum of the potential, and 
include all surrounding particles when computing profiles. 

The simple analytical description in (l)-(3) of the 
spherically averaged profiles of dark matter haloes has many 
applications. It will enable refined modelling of galaxy for- 
mation and dynamics, and improved calculations of gravi- 
tational lensing properties. A similarly detailed analysis of 
haloes formed in open models and models with a cosmolog- 
ical constant will be extremely interesting, and may lead to 
new diagnostics in which the halo density profiles can be 
used to discriminate between different cosmological models. 
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